REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-0188 

The  public  [ sporting  burden,  lor  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  Ihe  lime  lor  reviewing  instructions,  searching  existing  date  sources, 
gathering  and  maintaining  ihe  data  needed,  and  completing  and  reviewing  ihe  collection  of  Information.  Send  commems  regarding  this  burden  estimate  Of  any  other  aspect  ol  this  collection  of 
info* nation,  including  suggestions  lor  reducing  the  burden,  to  the  Dapanmeni  of  Defense.  Executive  Services  end  Communications  Directorate  j0?04  018BJ.  Respondents  should  be  aware 
that  notwithstanding  any  other  provision  of  law.  no  person  shall  be  sub[eti  to  any  penalty  lor  failing  to  comply  with  a  collection  of  information  if  rt  does  not  display  a  currently  valid  QMB 
control  n  urn  bet. 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ORGANIZATION. 

1.  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

18-07-2007  Journal  Article 

3.  DATES  COVERED  { From  -  To) 

4.  TITLE  AND  SUBTITLE 

A  Lagrangian  subgridscale  model  for  particle  transport  improvement  and 
application  in  the  Adriatic  Sea  using  the  Navy  Coastal  Ocean  Model 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

0602435N 

6.  AUTHOR(S) 

Angel ique  C  Haza;  Leonid  L  Piterbarg,  Paul  Martin,  Tamay  M.  Ozgdkmen, 
Annalisa  Griffa 

5d.  PROJECT  NUMBER 

5e,  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

73-6648-06-5 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Naval  Research  Laboratory 

Oceanography  Division 

Stcnnis  Space  Center,  MS  39529-5004 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

N  R  L/J  A/7 3 2 0-06-62 2 5 

9,  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Office  of  Naval  Research 

800  N.  Quincy  St, 

Arlington,  VA  22217-5660 

10.  SPONSOR/MONITOR  S  ACRONYM(S) 

ONR 

11.  SPONSOR/MONITOR’S  REPORT 

NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release,  distribution  is  unlimited. 


13.  SUPPLEMENTARY  NOTES 


14.  ABSTRACT 

An  accurate  estimation  of  Lagra  ngian  transport  in  the  ocean  is  important  for  a  number  of  practical  problems  such  as  dispersion  of  pollutants,  biological  species,  and 
sediments.  Forecasting  of  the  Lagmngian  pathways  necessarily  relies  on  the  accuracy  of  ocean  and  coastal  models.  However,  these  models  include  a  number  of  errors 
that  propagate  directly  from  the  Eulcrlan  velocity  Held  to  the  lagrangian  transport.  In  this  study,  so-called  Lagrangian  sub-grid-scale,  or  LSGS,  model  is  developed 
to  reduce  errors  projected  to  Lagrangian  transport  from  errors  arising  from  missing  phy  sics,  uncertainties  In  forcing  and  unresolved  scales  in  QGCMs.  The  LSGS 
method  acts  on  the  diagnostics  of  particle  transport  computed  from  coastal  or  ocean  models,  and  it  allow  s  to  minimize  the  discrepancy  between  the  statistical 
behavior  of  ihe  modeled  (synthetic)  and  real  (observed)  trajectories.  The  method  is  shown  to  work  well  using  both  a  so-called  Markov  velocity  field  model, 
representing  an  idealized  turbulent  flow  field,  and  in  the  context  of  the  Navy  Coastal  Ocean  Model  (NCOM)  configured  in  the  Adriatic  Sea  for  realistic,  high- 
resolution,  complex  ocean  flows.  The  simplicity  and  computational  efficiency  of  this  technique,  combined  with  applicability  to  ocean  models  at  a  wide  range  of 
resolutions,  appears  promising  in  light  of  the  challenge  of  capturing  exactly  the  oceanic  turbulent  fields,  which  is  critical  for  Lagrangian  dispersion. 


15.  SUBJECT  TERMS 

Lagrangian  transport;  pathways;  models;  trajectories 


16,  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

18.  NUMBER 
OF 

PAGES 

24 

19a,  NAME  OF  RESPONSIBLE  PERSON 

o.  REPORT 

b.  ABSTRACT 

C.  THIS  PAGE 

Paul  J.  Martin 

Unclassified 

Unclassified 

Unclassified 

UL 

19b.  TELEPHONE  NUMBER  (Include  area  coda) 
228-688-5447 

Standard  Form  298  (Rev.  8/98) 

Prescribed  by  ANSI  Std  Z39.18 


* 


Available  online  at  www.sciencedirectcom 


ScienceDirect 


Ocean 

Modelling 


ELSEVIER 


Ocean  Modelling  17  (2007)  68—91 


www.elsevier.coni/locate/occmod 


A  Lagrangian  subgridscale  model  for  particle  transport 
improvement  and  application  in  the  Adriatic  Sea 
using  the  Navy  Coastal  Ocean  Model 

Angelique  C.  Haza a,  Leonid  I.  Piterbarg  b,  Paul  Martin c, 

Tamay  M.  Ozgokmen a’*,  Annalisa  Griffaa,d 


*  RSMAS/MPO,  University  of  Miami,  Miami,  FLt  USA 
hCAMS.  University  of  Southern  California.  Los  Angeles.  CA ,  USA 
c  Naval  Research  Laboratory t  Stennis  Space  Center .  MS.  USA 
J  Consigtio  Nazionale  delie  Ricerche,  ISM  A  R.  LaSpezia,  hah 


Received  3  May  2006:  received  in  revised  form  27  October  2006;  accepted  30  October  2006 
Available  online  22  December  2006 


Abstract 


An  accurate  estimation  of  Lagrangian  transport  in  the  ocean  is  important  for  a  number  of  practical  problems  such  as 
dispersion  of  pollutants,  biological  species,  and  sediments.  Forecasting  of  the  Lagrangian  pathways  necessarily  relies  on 
the  accuracy  of  ocean  and  coastal  models.  However,  these  models  include  a  number  of  errors  that  propagate  directly  from 
the  Eulerian  velocity  field  to  the  Lagrangian  transport. 

In  this  study,  so-called  Lagrangian  sub-grid-scale,  or  LSGS,  model  is  developed  to  reduce  errors  projected  to  Lagrang- 
ian  transport  from  errors  arising  from  missing  physics,  uncertainties  in  forcing  and  unresolved  scales  in  OGCMs.  The 
LSGS  method  acts  on  the  diagnostics  of  particle  transport  computed  from  coastal  or  ocean  models,  and  it  allows  to 
minimize  the  discrepancy  between  the  statistical  behavior  of  the  modeled  (synthetic)  and  real  (observed)  trajectories. 
The  method  is  shown  to  work  well  using  both  a  so-called  Markov  velocity  field  model,  representing  an  idealized  turbulent 
flow  field,  and  in  the  context  of  the  Navy  Coastal  Ocean  Model  (NCOM)  configured  in  the  Adriatic  Sea  for  realistic,  high- 
resolution,  complex  ocean  flows. 

The  simplicity  and  computational  efficiency  of  this  technique,  combined  with  applicability  to  ocean  models  at  a  wide 
range  of  resolutions,  appears  promising  in  light  of  the  challenge  of  capturing  exactly  the  oceanic  turbulent  fields,  which 
is  critical  for  Lagrangian  dispersion, 

€>  2006  Elsevier  Ltd,  All  rights  reserved. 


1.  Introduction 

Ocean  general  circulation  models  (OCGMs)  arc  providing  increasingly  more  realistic  Eulerian  velocity 
fields  (e.g.,  Smith  et  ah,  2000;  Garraffo  et  al,  2001;  McClean  et  aL,  2002),  Subsequently,  there  has  been 
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a  development  of  operational  programs  that  rely  on  the  assimilation  of  various  observational  data  in  ocean 
circulation  models  to  forecast  the  ocean  state  (e.g.,  Pinardi  el  al.,  2003).  Realistic  ocean  and  coastal  model 
results  are  then  used  for  a  number  of  practical  applications,  such  as  ship  navigation  and  the  estimation  of 
the  dispersion  of  pollutants,  biological  species,  and  sediments.  At  the  basis  of  such  applications  is  the  trans¬ 
port  of  a  passive  tracer  field  and/or  particles.  It  is  this  particular  transport  problem  that  poses  significant 
computational  challenges  and  is  the  focus  of  the  present  investigation. 

Motivated  by  the  increase  in  the  number  of  drifters  and  floats  released  in  the  oceans  in  the  last  few  decades 
(e.g.,  Davis,  1991;  Owens,  1991;  Fralantoni,  2001:  Richardson,  2001;  Zhang  el  al.,  2001;  Lavender  el  al.,  2002; 
Bauer  el  al.,  2003;  Zhou  et  al.,  2002;  Reverdin  et  al.,  2003),  the  theoretical  treatment  of  the  Lagrangian  trans¬ 
port  problem  has  received  much  attention.  Two  general  theoretical  approaches  have  been  developed  to  under¬ 
stand  and  model  the  Lagrangian  transport.  The  first  is  a  statistical  approach,  which  is  motivated  by  the  work 
of  Thomson  (1986)  and  was  pioneered  for  oceanic  applications  by  Griffa  (1996).  In  this  approach,  the  two- 
dimensional  flow  field  is  decomposed  into  a  mean  and  eddy  velocity,  in  which  the  eddy  Lagrangian  velocity 
is  modeled  as  a  Markov  process  satisfying  the  so-called  Langevin  equation  or  random-flight  model.  This 
Lagrangian  stochastic  model  is  represented  by  a  first-order,  linear,  ordinary  differential  equation,  the  para¬ 
meters  of  which,  namely  the  velocity  variance  and  correlation  time,  can  be  estimated  using  drifter  data  sets. 
It  has  been  shown  by  Falco  et  al.  (2000)  that  such  models  can  be  useful  Tor  approximating  realistic  particle 
motion.  Other  versions  of  such  stochastic  models  have  been  developed,  namely  higher-order  models  (Berloff 
and  McWilliams,  2002),  multi-particle  models  (Piterbarg.  2001a,b),  and  models  that  take  into  account  the 
motion  due  to  trapping  in  mesoscale  vortices  (Veneziani  et  al.,  2005a, b).  Each  of  these  advances  requires  addi¬ 
tional  parameters  (for  instance,  the  spatial  decorrelation  scale  and  spin  parameter  in  the  latter  two  cases)  that 
need  to  be  estimated  from  oceanic  drifter  data.  Another  variation  is  to  assimilate  drifter  data  in  such  Lagrang¬ 
ian  models  to  help  predict  particle  motion  (Ozgokmen  et  al.,  2000,  2001;  Castellari  et  al.,  2001 ).  To  summarize, 
the  main  concept  in  these  studies  relics  on  a  stochastic  process  to  generate  simulated  trajectories  using  a 
model,  the  parameters  of  which  are  estimated  statistically  from  existing  data  sets.  One  of  the  fundamental 
problems  associated  with  this  approach  is  the  decomposition  into  a  mean  and  eddy  component,  which,  similar 
to  the  Reynolds  decomposition,  implies  that  the  eddy  component  represents  all  turbulent  scales.  Given  the 
high  Reynolds  number  of  oceanic  flows  exhibiting  coherent  structures  at  a  wide  range  of  scales  due  to  many 
interacting  processes,  the  underlying  assumption  of  Brownian  flow  in  such  models  is  idealistic. 

A  completely  different,  namely  purely  deterministic,  approach  has  been  developed  based  on  dynamical 
system  theory.  These  techniques  focus  on  identifying  coherent  structures  geometrically,  in  particular,  so-called 
hyperbolic  trajectories  in  the  Lagrangian  frame,  which  are  characterized  by  the  intersection  of  a  stable 
manifold  (along  which  fluid  particles  are  attracted  toward  the  hyperbolic  point)  and  an  unstable  manifold 
(along  which  fluid  particles  are  repelled  away  from  the  hyperbolic  point)  from  a  velocity  field  displaying 
complex  time  variability  (Wiggins,  2005).  One  of  the  issues  regarding  the  calculation  of  these  manifolds  is  that 
infinite  data  are  required  both  forward  and  backward  in  time,  which  is  clearly  impractical  for  oceanic  appli¬ 
cations,  Haller  and  Poje  (1998)  focused  on  the  transient  stagnation  points  and  finite-time  analogs  of  stable  and 
unstable  manifolds,  and  extended  the  lobe  analysis  to  the  treatment  of  finite-time  data.  This  method  was  then 
used  to  analyze  fluid  particle  pathways  in  an  eddy  resolving,  barotropic  model  (Poje  and  Haller,  1999).  Finite- 
time  geometric  techniques  have  been  successful  in  locating  the  boundaries  of  mesoscale  coherent  features  in 
the  Lagrangian  frame  (e.g.,  Coulliette  and  Wiggins,  2000;  Miller  et  al.,  2002).  It  has  also  been  shown  that  seed¬ 
ing  observations  in  rapidly  stretching  regions  of  the  flow  field  in  the  vicinity  of  hyperbolic  trajectories  leads  to 
improvement  of  the  reconstruction  of  Eulerian  fields  from  Lagrangian  data  (Poje  et  at.,  2002;  Toner  and 
Poje,  2004)  and  in  the  performance  of  a  Lagrangian  data  assimilation  scheme  (Moleard  et  al.,  2006),  However, 
only  persistent,  coherent  structures  can  be  captured  by  such  geometric  techniques  and  flow  regimes  with  no 
or  highly-transienl  coherent  structures  cannot  be  handled  effectively.  The  most  critical  aspect  of  this 
approach  is  the  total  reliance  of  the  Lagrangian  transport  map  on  the  accuracy  of  the  simulated  Eulerian 
velocity  fields. 

The  modeled  Eulerian  field  produced  by  the  OCGMs  cannot  be  entirely  accurate  because  of  three  funda¬ 
mental  reasons.  First,  the  model  forcing  via  wind  stress,  heat  flux,  precipitation,  evaporation,  and  boundary 
conditions  introduces  errors  due  to  space  and  time  sparseness  of  the  observational  data  sets  and  inaccuracies 
in  the  larger-scale  coupled  simulation.  Second,  ocean  and  coastal  models  typically  contain  missing  physics.  For 
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instance,  flow  interaction  with  capes  and  headlands,  which  are  ubiquitous  coastal  features,  can  introduce 
stratified  mixing  and  high  vertical  velocities  (e.g.,  Farmer  et  ah,  2002;  Pawlak  el  ah,  2003),  which  cannot 
be  handled  within  the  context  of  the  hydrostatic  approximation  generally  employed  in  such  models.  Surface 
mixed  layers  include  a  range  of  nonhydrostatic  stratified  mixing  processes  and  air-sea  interaction  problems 
(Kantha  and  Clayson,  2000)  that  pose  great  challenges  for  circulation  models.  Third,  given  that  circulation 
models  are  discretized  in  such  a  way  that  not  all  scales  of  motion  are  resolved  (typical  mesh  spacings  at  the 
present  time  are  100  km  in  climate  models,  10  km  in  general  circulation  models,  and  1  km  in  coastal  models), 
there  is  the  issue  of  the  effect  of  unresolved  or  filtered  scales .  For  instance,  most  OCGMs  use  Laplacian  or 
biharmonic  subgridscale  operators  for  momentum  and  active  tracers,  which  are  ultimately  based  on  the 
studies  by  Taylor  (1921)  and  Kolmogoroff  (1941)  of  laboratory-scale  flow  behavior.  There  are  no  oceanic 
observations  to  support  the  validity  of  such  a  generic  down-gradient  flux  assumption  or  to  accurately  estimate 
the  value  of  the  eddy  viscosity  and  diffusivity  as  a  function  of  the  mesh  spacing.  In  contrast,  the  development 
of  high-frequency  Doppler  radar  for  coastal  observations  revealed  submesoseale  (2-3  km  in  diameter)  eddies 
(Shay  et  al.,  2000),  the  parameterization  of  which  differs  significantly  from  that  based  on  eddy-viscosity 
models  (Caglar  el  al,,  2006).  Finally,  the  fractal  nature  of  the  Earth's  morphology  ( Weissel  el  al.,  1994)  implies 
unresolved  topography  at  even  the  smallest  scales. 

Thus,  it  is  clear  that  OCGMs  can  only  provide  an  approximate  Eulcrian  velocity  field  representing  the 
coastal  and  ocean  circulation.  Given  that  drifters  feel  forces  acting  on  the  scale  of  their  physical  size,  namely 
O  ( l  m)  (and  similar  considerations  for  other  passive  particles),  their  transport  with  the  flow  field  generated  by 
numerical  models  will  be  subject  to  a  range  of  errors  that  can  significantly  affect  their  paths  and  dispersion  (as 
implied  by  the  results  presented  in  Griffa  et  ak,  2004).  Since  such  circulation  models  are  computationally 
expensive  and  will  continue  to  contain  similar  errors  in  the  foreseeable  future,  how  can  their  Lagrangian  trans¬ 
port  characteristics  be  improved  in  a  cost-effective  and  realistic  manner? 

Here,  we  put  forth  and  explore  a  simple  method  as  a  first  step  toward  addressing  this  problem.  The  main 
underlying  idea  is  to  determine  and  minimize  the  discrepancy  between  the  statistical  behavior  of  the  modeled 
and  real  trajectories,  A  simple  model,  referred  to  as  a  Lagrangian  sub-grid-scale  (LSGS)  model  hereafter,  is 
then  developed  to  reduce  this  discrepancy.  The  method  is  based  on  the  use  of  stochastic  models,  but  differently 
from  previous  works  (e.g.  Griffa,  1996;  Berloff  and  McWilliams,  2002)  in  wrhich  the  focus  was  on  character¬ 
izing  the  action  of  the  entire  eddy  field,  here  the  final  goal  is  to  isolate  and  parameterize  only  that  part  of  the 
eddy  field  which  is  not  correctly  captured  and  reproduced  by  the  numerical  circulation  model.  The  approach  is 
inherently  statistical,  and  differs,  for  instance,  from  the  data  assimilation  approach,  in  which  Lagrangian 
observations  are  used  to  improve  prediction  of  a  specific  realization  of  velocity  or  trajectories  (e.g.,  Ozgdkmcn 
et  al.,  2000;  Molcard  ct  al.,  2003;  Taillandier  et  al,,  2006).  Here,  the  focus  is  on  particle  statistics  associated 
with  a  given  Eulcrian  velocity,  which  is  improved  using  information  from  observational  statistics.  The  LSGS 
model  is  capable  of  producing  ensemble  particle  trajectories  from  given  initial  conditions,  providing  informa¬ 
tion  on  tracer  dispersion  or,  equivalently,  on  the  maps  of  probability  of  finding  a  particle  at  a  given  point. 
Alternatively,  this  method  can  be  viewed  as  the  Lagrangian  counterpart  of  the  so-called  large-eddy-simulation 
method  (Sagaut,  2005),  in  which  the  main  concept  is  to  incorporate  the  effect  of  the  unresolved  turbulence  on 
the  resolved  turbulent  motion  by  adding  SGS  stresses  to  the  Eulcrian  momentum  equations.  These  SGS  stres¬ 
ses  are  tuned  to  be  consistent  with  physical  insight  (e.g.,  the  Kolmogoroff  turbulent  energy  spectrum)  or  with 
results  from  fully-resolved  simulations.  In  our  case,  the  LSGS  model  represents  the  terms  added  to  the 
Lagrangian  advection  model  such  that  particle  advection  is  modified  to  be  consistent  with  the  behavior  of 
observed  trajectories  in  the  statistical  sense. 

The  method  is  tested  in  three  stages.  First,  a  Lagrangian  stochastic  model  in  the  context  of  the  Langcvm 
equation,  a  random-flight  model,  is  employed  to  model  the  Lagrangian  eddy  field  to  confirm  that  the  LSGS 
model  acts  to  restore  the  statistical  properties  of  the  observed/reference  particle  motion.  The  random-flight 
model  generates  particle  trajectories  based  on  input  statistical  parameters,  such  as  the  turbulent  velocity  fluc¬ 
tuation  and  correlation  time  scale,  without  the  need  of  an  underlying  Eulcrian  velocity  field.  Thus,  it  is  a  first 
ideal  testbed.  However,  the  correction  by  the  LSGS  model  in  the  context  of  the  random-flight  model  does  not 
take  into  account  the  effect  of  the  local  model  fields.  Thus,  a  so-called  Markov  velocity-field  model,  represent¬ 
ing  the  behavior  of  an  idealized,  two-dimensional,  turbulent  velocity  field,  is  employed  to  test  a  space-depen¬ 
dent  version  of  the  method.  Finally,  the  LSGS  model  is  applied  to  rectify  the  particle  transport  based  on 
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hourly  output  from  a  realistic,  high-resolution  ( 1-krn)  ocean  model,  the  Navy  Coastal  Ocean  Model  (NCOM), 
configured  in  the  Adriatic  Sea.  The  LSGS  technique  offers  two  significant  advantages.  The  first  is  the  funda¬ 
mental  issue  that  it  can  be  applied  to  numerical  models  at  all  resolutions,  since  the  principle  of  the  LSGS 
model  is  to  correct  the  errors  in  the  Lagrangian  transport  due  to  the  effect  of  unresolved  scales  of  motion, 
missing  physics,  and  forcing  errors.  The  second  is  that  the  LSGS  model  can  be  used  on  the  transport  diagnos¬ 
tics  from  post-processed  (low-pass  filtered)  coastal  or  ocean  model  outputs.  In  other  words,  existing  Eulerian 
output  can  be  employed  in  conjunction  with  this  technique  in  order  to  obtain  a  more  realistic  particle  trans¬ 
port,  without  the  need  to  rerun  the  circulation  model.  The  main  requirement  (and  possible  drawback)  is  that 
an  observational  drifter  data  set  of  adequate  coverage  (space/time  density)  is  needed  to  improve  the  modeled 
Lagrangian  transport  in  the  area  of  interest. 

The  paper  is  organized  as  follows:  The  theory  leading  to  the  LSGS  is  provided  in  Section  2.  The  proof  of 
the  LGSG  model  formula  and  details  of  the  Markov  velocity  field  model  are  given  in  Appendix  A.  NCOM  is 
described  in  Section  3.  The  results  from  both  the  Markov  velocity  field  model  and  NCOM  are  presented  in 
Section  4.  The  conclusions  and  directions  for  future  work  are  discussed  in  Section  5, 

2.  Theory  of  the  Lagrangian  subgridscale  model 

Lagrangian  motion  in  a  random  or  even  complex  deterministic  velocity  field  can  be  described  as  a  stochas¬ 
tic  process.  Thus,  its  full  description  is  given  by  a  joint  probability  distribution  of  velocities  and  positions  of 
any  number  of  particles  for  any  set  of  time  moments.  It  is  not  realistic  to  estimate  or  derive  theoretically  this 
distribution  in  any  region  of  the  world  ocean  for  now  or  in  the  future.  To  this  end,  only  the  statistical  char¬ 
acteristics  of  one-particle  motion  have  been  extensively  studied  (e.g.,  Veneziani  et  al.,  2004;  and  references 
cited  therein)  and  some  efforts  have  been  carried  out  toward  two-particle  motion  (e.g.,  LaCasce  and  Bower, 
2000;  LaCasce  and  Ohlmann,  2003).  Thus,  in  this  paper  we  focus  on  the  statistics  of  single-particle  motion 
characterized  by  a  joint  pdf  p(/,v,  r)  of  the  velocity  v  and  position  r  at  moment  t.  Recall  that  the  one-particle 
motion  pdf  completely  determines  the  mean  concentration  c(r)  of  a  passive  scalar,  while  the  two-particle  sta¬ 
tistics  allow  finding  the  spatial  covariance  tensor  of  the  concentration  field.  In  particular,  for  an  incompress¬ 
ible,  inviscid  flow  and  a  delta-function  source,  c(r)  is  obtained  by  integration  of  the  pdf  over  v. 

In  turn,  an  OGCM  produces  an  ensemble  of  Lagrangian  trajectories  characterized  by  the  model  pdf 
pm(/,v,r).  For  numerous  reasons  indicated  in  Section  1,  pm(r,v,r)  turns  out  to  be  quite  different  from 
p(t,  v,  r).  Therefore  a  typical  OGCM  is  not  able  to  reproduce  dispersion  of  a  passive  scalar  with  high  accuracy. 
A  natural  question  is  whether  a  model  trajectory  ensemble  can  be  corrected  in  a  way  to  approximately  retrieve 
the  real  pdf  p(t,  v,r).  The  most  obvious  correction  would  be  the  addition  of  another  stochastic  process  to  the 
model  velocity  to  compensate  for  the  effect  of  missing  physics  and  scales.  Thus,  the  simplest  mathematical 
formulation  of  a  LSGS  model  is  as  follows. 

Problem  1.  Let  v(r),  vm(r)  be  the  real  and  model  sample  Lagrangian  velocities ,  respectively,  of  a  single  particle  in 
the  lime  interval  [0,7],  Find  a  random  vector  process,  the  SGS  or  missing  velocity  component,  >}(t),  such  that  the 
corrected  position  and  velocity  defined  by 

dtc{t)/dt  =  vc(r),  vc(/)  =  vm(r)  +  i/(/}  (1) 

have  the  same  pdf  p(t.v,r)  as  a  real  particle  provided  with  the  same  initial  condition,  for  each  time  moment 
t  €  [0,  71. 

Thus,  by  inserting  the  missing  component,  one  can  improve  the  prediction  of  tracer  spreading. 

In  a  significant  part  of  the  world  ocean,  individual  Lagrangian  trajectories  on  horizontal  or  isopycnal 
surfaces  can  be  described,  at  least  to  first  approximation,  by  the  first-order  Markov  or  random-flight  model 
(see  again  Veneziani  et  al,,  2004;  Faleo  el  al.,  2000) 

dr  =  (V (r.  r)  +  V)  dt ,  dv'  =  -AV  dr  +  A  d w,  (2) 

where  V(r,r(/))  is  the  deterministic  drift,  v'  is  the  fluctuation  velocity  with  zero  mean  driven  by  standard  2D 
Brownian  motion  w.  The  dissipation  and  dispersion  matrices  are  given  by 
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(3) 


where  o,„  ov  are  the  rms  of  the  zonal  and  meridional  velocity  fluctuations,  respectively,  iu,  t„  are  the  corre¬ 
sponding  Lagrangian  correlation  times,  and  Q  is  the  spin. 

The  model  was  first  developed  for  Lagrangian  motion  in  the  atmosphere  by  Thomson  (1986),  and  the  con¬ 
cept  of  spin  was  introduced  by  Borgas  et  al,  (1997)  and  Reynolds  (2002,  2003).  A  remarkable  fact  is  that  in  the 
framework  of  (2),  the  aforementioned  pdf  p(t.  v,  r)  is  completely  determined  by  the  initial  conditions,  the  drift 
V(f,r),  the  turbulent  velocity  variance  a,  the  Lagrangian  correlation  timer,  and  the  spin  Q.  More  exactly,  this 
pdf  can  be  found  as  a  solution  of  the  corresponding  Fokker-Planck  equation  given  the  listed  parameters. 

The  basic  assumption  behind  our  approach  is  the  following: 


Assumption  I.  Both  the  real  and  model  single-particle  trajectories  are  well  approximated  by  model  (2,3)  with 
the  same  drift  V  and  different  fluctuation  parameters  a,  r,  Q. 


Under  this  assumption,  which  in  particular  implies  that  the  model  is  capable  of  exactly  capturing  the  drift. 
Problem  l  can  be  reformulated  as  follows: 

Problem  2.  Let  v'f  (/)  and  (/)  be  the  real  and  model  Lagrangian  jluctuation  velocities  of  a  single  particle , 
respectively,  and  described  by  the  same  Langevin  equation  in  (2)  with  different  parameters.  Find  a  random  vector 
process  (missing  component)  >;(t)  such  that  the  corrected  velocity  fluctuation 

<(0  =  +  9(0  ‘  (4) 

is  covered  by  the  Langevin  equation  with  real  parameters. 

In  other  words,  we  want  the  corrected  velocity  and  position  to  satisfy  (2)  with  real  parameters,  but  of  course 
with  another  (say  independent)  Brownian  forcing  since  the  goal  is  to  reproduce  just  the  real  statistics  rather 
than  an  individual  trajectory  itself  This  problem  formulation  requires  knowledge  of  the  real  fluctuation 
parameters.  In  addition,  it  is  natural  to  require  the  missing  component  to  be  completely  determined  only 
by  the  model  trajectory  (plus,  of  course,  the  statistics  of  the  real  trajectories). 

A  solution  of  the  formulated  problem  is  given  by  the  following  statement*  Introduce 

Cr(x)  =  x I  +  Ar>  Cm(x)  =  xl  +  Am>  Q(x)  =  AtCt(*)A;\  P(x)  =  ArCm{x)A„'  -  Q(x), 

where  matrices  Ar,  At  and  Am,  Am  are  defined  by  (3)  for  real  and  model  parameters,  respectively,  /  is  the  2  x  2 
unit  matrix,  and  x  is  a  dummy  variable. 

Proposition  1.  Let  the  2D  process  i)(t)  be  a  solution  of  the  following  stochastic  differential  equation 

e{i)"=p{^  (5) 


Then,  the  corrected  velocity  defined  by  (4)  satisfies  the  equation  for  the  real  velocity  in  (2),  thereby  yielding  a 
solution  of  Problem  2. 


The  proof  is  given  in  Appendix  A.  1 .  If  both  spins  arc  zero,  Qr  =  f2m  =  0,  then  (5)  breaks  down  into  two 
separate  identical  equations  for  each  component.  The  equation  for  the  zonal  direction  becomes: 


d»f 

d/ 


lit 


+  bu'm  +  «|, 


(6) 


where 


^  r  y/^m  j 


b  = 


7 

Tr 


C 


(?) 


Now  subscripts  r,  m  denote  real  and  model,  respectively.  In  particular,  the  variance  (energy)  of  the  missing 
component  >j  in  the  stationary  case  is  given  by 
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Tr  +  Tm 

which  is  zero  if  the  model  and  real  parameters  coincide. 

The  correlation  coefficient  rmc  between  the  modeled  and  corrected  velocity  also  can  be  found  from  (6),  (7) 
(see  Appendix  A. 2),  It  is  interesting  that  it  depends  only  on  the  ratio  y  —  tni/rr  of  the  model  and  real  Lagrang- 
ian  time  scales  and  has  a  simple  analytical  expression,  namely 


In  particular,  for  1/4  ^  y  ^  4  we  have  rmc  ^  0,8.  That  high  correlation  is  well  illustrated  by  the  experiments 
following  below. 

To  ensure  stationarity  of  the  LSGS  component  >/(*),  the  initial  condition  jj(0)  should  be  taken  as  a  normal 
random  value  with  zero  mean  and  with  variance  given  by  (8),  For  long  enough  realizations,  i/(f)  becomes 
stationary  for  /  tr  for  arbitrary  initial  conditions. 

To  illustrate  how  the  suggested  procedure  of  adding  the  LSGS  component  performs,  we  first  experiment 
with  a  simplified  version  of  the  model  (2),  in  which  the  mean  drift  is  zero  and  the  variances  of  the  zonal 
and  meridional  components  arc  equal,  as  well  as  the  correlation  times,  for  both  the  real  and  model  flows. 
The  scheme  of  the  simulations  is  as  follows.  Given  a  set  of  the  real  aTt  Qr  and  model  am,  rmT  Qm  parameters, 
the  real  and  model  velocities  and  corresponding  trajectories  are  generated  via  (2),  Then,  by  solving  (6)  or  (5) 
(if  the  spin  is  not  zero),  the  LSGS  component  is  computed  and  inserted  by  (I),  Finally,  the  corrected  and  real 
trajectories  (velocities)  are  compared  visually  and  statistically. 

In  Fig.  1,  an  example  of  the  model  real  and  corrected  velocities  (  u-componcnt)  arc  shown  for  an  integration 
time  period  of  T  =  100  days  and  parameter  values  aT  =  20 cm/s,  rr  =  6days,  Qv  =  0,  <7m=10cm/s, 
rm  —  1.5  days,  and  Qm  =  0,  Since  the  spin  is  equal  to  zero,  we  used  Eq,  (6),  It  can  be  seen  that  the  corrected 
velocity  has  energy  values  and  fluctuation  scales  closer  to  the  real  ones  with  respect  to  the  model  velocity, 
while  it  is  well  correlated  with  the  latter.  The  parameters  estimated  from  the  corrected  velocity. 


Real(top).  Modef(mid),  Correct ed( bottom )P  am=10h  if=4imi  or=2om 


days 

Fig.  I  Time  series  of  model  (red),  real  (blue),  and  corrected  (black)  velocities.  The  lime  series  of  the  real  (corrected)  velocities  is  shilled  by 
50  cm  s" 1  upwards  (downwards)  lo  belter  demonstrate  their  differences.  In  particular  note  that  the  LSGS  model  acts  to  enhance  the 
amplitude  of  the  model  velocity  fluctuations  to  be  more  consistent  with  the  real  case.  The  means  of  all  three  time  series  are  zero, 
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Trajectories;  real(biue),  model(redJf  corrected  (black) 


Fig.  2.  Model  (red),  real  (blue),  and  corrected  (black)  trajec lories. 


<rc  —  20  8  cm/s  and  xc  =  5*  16  days,  are  dose  to  the  real  values.  The  corresponding  model  trajectory  is  much 
shorter  than  the  real  one  due  to  the  deficit  of  energy  while  the  corrected  one  is  approximately  of  the  same 
length  as  the  real  trajectory  (Fig.  2). 

Another  experiment  with  smaller  velocities  and  non-zero  real  spin  with  trt  =  1  cm/s,  tr  -  3.5  days, 
£2r  =  2n/20,  =  0.25  cm/s,  rm  —  10  days,  and  Qm  =  0  shows  that  the  LSGS  algorithm  is  able  to  reproduce 

loops  observed  in  the  real  trajectory  even  though  there  are  no  loops  in  the  model  trajectory  (Fig.  3). 

Notice  that  so  far,  there  is  no  explicit  spatial  structure  of  the  motion  involved  in  our  approach.  The  LSGS 
procedure  given  by  (6)  or  (5)  is,  in  fact,  a  posterior  procedure,  since  a  Lagrangian  trajectory  is  first  computed 
over  the  whole  time  interval  [0,7]  and  then  Eq.  (6)  or  (5)  is  solved  in  the  same  interval.  Thus,  no  adjustment  to 
the  factual  position  of  the  particle  is  suggested,  which  can  be  misleading  if  the  particle  motion  in  a  velocity  field 
is  considered.  As  an  alternative,  we  present  the  so  called  space-dependent  LSGS  model  as  follows. 

First,  we  introduce  the  model  Eulerian  velocity  field  (say  zonal  component)  and  assume  a  constant  drift 

=  U  +  i4(/,r). 


The  assumption  is  due  to  difficulties  arising  from  the  separation  of  Lagrangian  motion  into  deterministic  drift 
and  fluctuations  for  a  non-constant  mean  flow,  say  U(f,  r).  More  exactly,  the  problem  is  that  in  this  case  the 
drift  V(f.  r)  introduced  in  (2)  is  not  equal  to  U(/,  r)  at  all,  and  moreover,  cannot  be  expressed  explicitly  in  terms 
of  U  and  the  statistics  of  u'.  Then  rewrite  the  Eqs.  (6)  and  (1)  for  the  corrected  velocity  and  zonal  component 
of  position  in  a  slightly  different  form: 


d>?(0  _  di4(/.rm(/)) 
d/  1  d  / 


dr 


(9) 


The  suggested  space-coordinate  dependent  procedure  is  obtained  by  replacing  rm(f)  in  (9)  with  rc(r),  i.c,: 


d>/(0 

dr 


dM^,(/,rc(/)) 

dr 


+  bum{t,  rc(/))  +  «/(/), 


d.vc(r) 

dr 


V  +  u'm{t,Tc(t))  + 


(10) 


This  formula,  together  with  an  obvious  replica  for  the  meridional  component,  is  a  basis  for  implementing  the 
LSGS  model  in  a  space-dependent  model.  In  the  following,  wre  will  test  it  using  both  a  synthetic  random  velo¬ 
city  field  and  NCOM,  While  proof  for  the  random-flight  model  is  given  in  Appendix  A.l,  we  cannot  present  a 
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Trajectories:  reai(biue)r  model(red),  corrected  (black) 


Pig,  3,  Model  I  red),  real  (blue),  and  corrected  (black)  trajectories,  in  the  case  where  missing  spin  is  added  by  the  LSGS  model. 

proof  for  the  ease  of  a  random  Markov  velocity  field  model  that  the  statistics  of  the  corrected  trajectories  (via 
Eq.  (10)  in  this  ease)  are  identical  to  those  of  the  real  trajectories-  It  is  probably  a  hard  mathematical  problem 
since  an  analysis  of  the  space-dependent  procedure  requires  taking  into  consideration  the  spatial  structure  of 
the  velocity  field-  For  this  reason,  the  procedure  has  been  thoroughly  tested  with  numerical  simulations.  As 
will  be  seen  later,  the  conjecture  will  turn  out  to  be  well  supported  by  the  simulations. 

3.  Navy  Coastal  Ocean  Model 

The  ocean  model  used  for  this  study  is  the  Navy  Coastal  Ocean  Model  (NCOM)  as  described  in  Martin 
(2000),  with  some  improvements  as  described  in  Morey  et  al,  (2003)  and  Barron  et  at.  (2006).  This  is  a  hydro¬ 
static  model,  which  is  similar  in  its  physics  and  numerics  to  the  Princeton  Ocean  Model  (POM)  (Blumberg  and 
Mellon  1987),  but  uses  an  implicit  treatment  of  the  free  surface  and  a  hybrid  vertical  grid  with  sigma  coordi¬ 
nates  in  the  upper  layers  and  (optionally)  level  coordinates  below  a  user-specified  depth. 

The  model  equations  include  a  source  term  that  can  be  used  for  river  inflows.  A  third-order  upwind  method 
(Holland  et  al.,  1998)  was  used  for  advection.  Vertical  mixing  was  computed  using  the  Meltor-Yamada  Level 
2  scheme  (Mellor  and  Yamada,  1974).  The  equation  of  state  used  was  that  of  Mellor  (1991). 

The  ocean  model  domain  consists  of  the  entire  Adriatic  Sea,  a  sub-basin  of  the  Mediterranean  Sea,  and  it 
includes  the  Strait  of  Otranto  and  a  small  part  of  the  northern  Ionian  Sea.  The  horizontal  grid  resolution  is 
1019.5  m.  The  vertical  grid  consists  of  32  total  layers,  with  22  sigma  layers  used  from  the  surface  down  to  a 
depth  of  291  m  and  level  coordinates  used  below  291  m.  Hence,  the  grid  is  like  a  regular  sigma  coordinate  grid 
in  water  shallower  than  291  m  and  similar  to  a  level  grid  in  deeper  water.  The  vertical  grid  is  uniformly 
stretched  from  the  surface  downward  with  a  maximum  thickness  of  the  upper  layer  of  2  m  and  a  maximum 
depth  of  1262  m. 

Initial  conditions  and  daily  boundary  conditions  (BC)  were  taken  from  a  hindcast  of  a  global  model 
(Barron  el  al..  2004),  The  numerical  treatment  of  the  BC  includes  the  Flathcr  radiation  condition  (Flalher 
and  Proctor,  1983)  for  the  surface  elevation  and  depth-averaged  normal  velocity,  Orlanski  radiation  condi¬ 
tions  (Orlanski,  1976)  for  the  tangential  velocities  and  scalar  fields,  and  a  relaxation  to  the  temperature 
and  salinity  from  the  global  model  near  the  open  boundary. 
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Tidal  forcing  was  provided  using  tidal  elevation  and  depth-averaged  normal  and  tangential  velocities  at  the 
open  boundaries  from  the  Oregon  State  University  (OSU)  tidal  data  bases,  which  are  derived  from  satellite 
altimetry  data  (Egbert  and  Erofeeva,  2003).  Data  from  the  OSU  Mediterranean  tidal  data  base  were  used 
for  the  Kl,  Ol,  M2,  and  S2  constituents  and  data  from  the  OSU  global  data  base  were  used  for  PI,  Ql, 
K2,  and  N2.  Tidal  potential  forcing  for  these  eight  constituents  was  used  in  the  interior  of  the  model  domain. 

Atmospheric  forcing  was  obtained  from  the  Coupled  Ocean/ Atmosphere  Mcsoscale  Prediction  System 
(CO AMPS)  (Hodur,  1997).  The  COAMPS  setup  for  the  Adriatic  consists  of  a  triply-nested  grid  with  resolu¬ 
tions  of  36,  12,  and  4  km  (Pullen  et  aL  2003).  The  outer  grid  of  this  nested  grid  system  covers  most  of  Europe 
and  the  Mediterranean  and  the  inner  4-kni  grid  covers  the  entire  Adriatic  and  part  of  the  Tyrhennian  Sea. 
COAMPS  itself  is  nested  within  the  Navy  Operational  Global  Atmospheric  Prediction  System  (NOGAPS) 
(Rosmond  et  ah,  2002). 

River  and  runoff  inflows  for  the  Adriatic  were  taken  from  the  monthly  climatological  data  base  of  Raicich 
(1994).  This  data  base  includes  discharges  for  about  39  rivers  and  runoff  inflows  along  a  number  of  sections  of 
the  Adriatic  coastline.  Raiciclrs  monthly  climate  values  were  used  for  all  the  inflows,  except  that  daily 
observed  discharge  values  were  used  for  the  Po  River  (Rich  Signed,  personal  communication). 

4.  Results 

4 .  L  General  experimental  strategy 

The  main  experimental  strategy  is  as  follows:  first,  the  space-dependent  LSGS  model  (10)  is  tested  within 
the  framework  of  a  purely  random  Eulerian  velocity  field  that  is  Markovian  in  time.  In  the  simplest  ease  of 
zero  mean  flow,  this  model  is  characterized  by  only  three  parameters:  the  velocity  variance,  the  correlation 
time,  and  the  spatial  correlation  radius.  For  this  velocity  field,  individual  Lagrangian  trajectories  are  described 
by  the  random-flight  model  (2)  with  a  high  degree  of  accuracy  (but  not  exactly).  Since  the  theory  is  based 
precisely  on  the  Lagrangian  equations  of  motion  given  by  (2),  it  is  expected  that  the  LSGS  would  act  to  make 
an  accurate  correction  toward  the  real  trajectory  properties. 

Then,  a  much  more  demanding  test  is  conducted  in  the  context  of  NCOM  in  the  Adriatic  Sea,  which  is 
subject  to  highly-variable  surface  wind  forcing,  contains  tides  and  surface  density  gradients.  In  addition, 
the  three-dimensional  geometry  of  the  domain  has  a  first-order  effect  on  the  How  field,  such  that  there  are 
a  rich  variety  of  co-existing  turbulent  flow  regimes.  Since,  generally  speaking,  Eq.  (2)  is  not  an  exact  model 
for  particles  transported  by  such  a  complex  flow  field,  the  theory  will  not  be  strictly  valid.  The  main  idea  is 
to  explore  how  the  LSGS  method  would  work  under  such  demanding  circumstances. 

The  main  experimental  matrix  common  to  both  the  Markov  velocity  field  model  and  NCOM  consists  of 
experiments  where  the  velocity  fluctuation  of  the  real  drifters  is  varied  in  the  range  of  —  L  =  l, 
^  —  2,  and  the  correlation  time  is  such  that  f  -  =  f-  =  1,^  =  4  (Table  1).  To  clearly  identify  the  changes 
due  to  the  LSGS  correction,  M  —  121  trajectories  are  released  in  each  experiment  with  different  initial  condi¬ 
tions  i/(0)  for  the  LSGS  component,  and  the  dispersion  is  estimated  from 

/ 1  m  \1/2 

(^E6m(0-rm(0))2J  •  (11) 

Table  1 


The  main  experimental  matrix 
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II  - 1 

Exp-I 
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Il  =  4 

Exp-6 

Exp-7 
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in 


A  C  Haza  et  ai  !  Ocean  Modelling  17  (2007)  68-91 


11 


Since  p{i)  ^  Gyfri^  changing  a  and  r  by  factors  of  two  and  four,  respectively,  should  lead  to  comparable 
impacts  on  the  trajectories.  A  large  factor  of  two,  rather  than  a  subtle  difference,  is  chosen  to  better  illustrate 
the  changes  due  to  LSGS  correction.  Thus,  eight  sets  of  LSGS  experiments  are  conducted  (the  case  with  f-  —  1 
and  =  1  requires  no  correction). 

4.2,  Results  from  the  Markov  velocity  field  model 

In  this  section,  the  results  are  reported  for  a  homogeneous,  stationary,  stochastic  incompressible  velocity 
field  with  zero  mean  and  with  the  covariance  of  the  stream-function  given  by 

Bv(t,r)  =  exp{  — |(|/te  —  r2  /R2},  (12) 

where  R  is  the  velocity  correlation  radius.  The  values  of  the  base  parameters,  (7m  =  lOcms  1  and 
te  =  1.5  days,  are  in  the  general  range  of  those  estimated  from  drifter  trajectories  in  the  Adriatic  Sea  (e.g,, 
Falco  et  at.T  2000;  Poulain,  2001),  However,  it  is  not  our  objective  to  try  to  produce  simulated  trajectories  that 
are  identical  to  those  in  the  Adriatic  Sea.  The  numerical  value  R  =  20  km  used  in  the  simulations  is  a  generic 
first-mode  radius  of  deformation.  The  mean  flow  field  is  assumed  to  be  zero  to  prevent  coherent  structures 
such  as  jets  and  recirculations  from  influencing  the  Lagrangian  statistics. 

The  field  (12)  is  Markovian  in  time.  However,  Lagrangian  trajectories  generated  by  it  arc  not  covered  by  (2) 
exactly.  Instead,  Eq.  (2)  yields  just  a  very  accurate  description  of  single-particle  motion  because  the  Kubo  non¬ 
linearity  parameter  K  =  umi E//?  =  0.648  is  relatively  small.  Numerous  experiments  with  different  R  showed 
that  for  K  <  1,  not  only  is  the  Lagrangian  velocity  well  approximated  by  the  Langevin  equation,  but  its  cor¬ 
relation  time  tm  is  very  close  to  the  Eulerian  correlation  time  rE.  For  large  K ,  this  is  not  the  case  and  the  depen¬ 
dence  of  rm  on  R  becomes  strong. 

In  the  first  series  of  experiments,  we  implemented  both  posterior  (6)  and  space-dependent  (10)  LSGS 
procedures  with  a  two-fold  goal:  first,  to  compare  the  two  methods,  and  second,  to  check  how  well  the  real 
statistics  are  reproduced  by  both  for  the  wide  range  of  parameters  given  in  Table  1. 

Eq.  (10)  were  solved  by  using  the  spectral  decomposition  of  the  velocity  field  by  making  use  of  the  Markov 
property  (see  Appendix  A, 2).  The  integrations  arc  carried  out  with  a  time  step  of  0.1  days  for  T  =  15  days, 
which  is  an  order  of  magnitude  larger  than  the  correlation  time  scale.  The  first  important  conclusion  from 
the  8  simulation  experiments  is  that  the  estimates  of  the  key  parameters  i  and  a  are  almost  the  same  for 
the  posterior  (6)  and  space-dependent  (10)  versions.  We  do  not  show  the  exact  numbers,  but  rather  state  that 
in  most  cases,  the  difference  between  the  parameters  rrc,  tc  of  the  corrected  velocity  and  the  corresponding  real 
parameters  does  not  exceed  5%.  We  point  out  again  that  this  result  is  due  to  small  Kubo  numbers  (or  relatively 
large  space-correlation  radii)  and,  in  the  opposite  case,  the  divergence  between  the  two  is  significant.  Hence, 
for  homogeneous  environments  and  small  Kubo  numbers,  the  space-dependent  LSGS  model  has  the  same 
property  as  the  posterior  procedure:  it  accurately  reproduces  the  real  statistics. 

Let  us  demonstrate  the  changes  in  dispersion  caused  by  the  LSGS  correction.  First,  spagetti  plots  are  shown 
in  Fig.  4  for  all  8  experiments.  A  total  of  121  particles  are  released  in  each  experiment  from  within  a 
10  km  x  10  km  area  with  a  spacing  l  km,  from  which  40  trajectories  are  depicted.  Besides  the  initial  positions, 
the  corrected  trajectories  differ  among  themselves  by  the  initial  conditions  for  the  LSGS  components  ij(0). 
Specifically,  the  ^(Q)’s  are  sampled  from  the  normal  distribution  with  zero  mean  and  variance  given  by  (8), 
All  the  trajectories  in  each  spagetti  diagram  are  obtained  from  the  same  velocity  field  realization  (12).  The  dia¬ 
grams  clearly  demonstrate  the  changes  in  the  corrected  trajectories  caused  by  the  changes  in  the  real  parameters. 

The  quantitative  comparison  between  the  targeted  and  achieved  statistical  parameters  (cr,  t)  in  Exp-1  to 
Exp-8  is  demonstrated  in  Fig.  5.  This  scatter  plot  shows  that  the  parameters  of  the  LSGS  experiments 
and  rc  (for  definiteness:  we  used  the  estimates  coming  from  the  space-dependent  version  (10))  are  successfully 
modified  to  change  from  (<rro,  Tni)  toward  (oy,  tr).  A  satisfactory  agreement  is  obtained  in  all  the  experiments* 
even  though  in  some  cases,  as  in  Exp.  1,  the  relative  error  can  reach  30-40%.  We  suggest  that  a  significant 
difference  between  <jc  and  <rr  in  this  case  could  be  caused  by  the  combination  of  two  factors:  first,  the  sample 
variability  of  the  statistical  estimate  of  the  corrected  velocity  and,  second,  imperfect  performance  of  the  LSGS 
correction.  Even  though  the  total  sample  size  N  —  ill  x  150  =  18, 150  looks  impressive,  the  equivalent  num¬ 
ber  of  independent  observations,  Ac,  is  much  smaller.  Indeed,  for  each  trajectory,  ~  15/ir  is  only  40  and 
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Fig,  4,  Fiftccn-day  long  trajectories  of  121  particles  for  Exp- 1  to  Exp-8  and  Set- A  generated  by  the  random  field  model  using  the  space- 
dependent  version  of  the  LSGS  model  Red  and  blue  points  mark  the  initial  and  final  positions*  respectively.  The  model  run  (without 
LSGS)  is  depicted  in  the  middle  panel. 


the  trajectories  in  each  duster  arc  strongly  dependent  since  they  start  dose  to  one  another,  which  probably 
makes  total  Ne  two  orders  of  magnitude  less  than  N.  As  for  the  correction  method  itself,  recall  that  the  theory 
presented  in  Section  2  is  exact  for  the  posterior  version  only  and  it  is  essentially  based  on  continuous  time 
consideration,  stationarity,  and  Gaussianity.  In  experiments,  the  first  condition  cannot  be  fullfillcd  while 
the  others  hold  only  approximately.  It  is  not  surprising  that  superposition  of  all  the  listed  factors  can  some¬ 
times  affect  results,  even  though  the  effect  of  each  one  of  them  seems  negligible. 

The  plots  of  p(t)  defined  in  (11)  and  those  obtained  from  the  basic  model  (2)  with  the  real  parameters 

P,U)  =  2<rr(Tr{f  -  Tr(l  -  exp(-f/Tr))))l/;! 

are  shown  in  Fig.  6  for  a  subset  of  experiments,  namely  for  Exps.  3*  4,  and  7,  Each  of  these  experiments  is 
taken  from  a  different  column  and  row  of  Table  1  to  better  represent  the  whole  variety  of  relations  between 
the  model  and  real  parameters.  The  model  curve  pm(f)  of  the  same  analytical  form  as  the  model  parameters 
is  also  plotted.  For  Exps.  3  and  4,  the  agreement  between  the  curves  is  very  good  in  both  the  initial  and 
final  stages.  For  Exp.  7  that  difference  is  significant.  Even  though  the  corrected  and  real  parameters  in  this 


AC  Haza  el  alt  Ocean  Modelling  17  (2007)  68-91 


79 


Fig.  5.  Scalier  plot  of  targeted  and  achieved  values  of  the  statistica  l  parameters  from  real  and  corrected  parlide  trajectories  in  experiments 
shown  in  Fig.  4.  The  numbers  in  black  denote  experiments  (Table  1 J  at  locations  white  numbers  in  green  mark  ,  £).  The  model 

parameters  are  marked  by  the  symbol  "m", 
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Fig,  6.  Time  evolutions  of  pi/)  from  the  random  field  model  (solid  lines)  and  pr(/)  (dashed  lines)  for  selected  experiments. 
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experiment  are  closer  than  those  in  Exps.  3  and  4,  the  divergence  at  the  15th  day  is  essentially  larger  because  of 
the  different  behavior  during  the  initial  stage.  In  Exps.  6  and  S,  the  gap  between  the  corrected  and  real  vari¬ 
ances  is  larger  than  those  in  Exps.  3,  4,  and  7  and  the  relative  differences  between  p(15)  and  pr(15)  arc  20%  and 
27%,  respectively.  In  summary,  the  method  not  only  gives  rise  to  the  real  parameters,  but  is  also  able  to  repro¬ 
duce  the  whole  time  evolution  of  the  dispersion  curve  satisfactorily. 

43,  Application  to  NCOM  in  the  Adriatic  Sea 

The  mean  surface  flow  field  from  the  NCOM  simulation  is  depicted  in  Fig.  7a.  The  reader  is  referred  to 
Poulain  (1999,  2001),  Cushman-Roisin  et  al.  (2001),  Falco  et  al.  (2000)  and  Maurizi  et  al,  (2004)  for  a  com¬ 
prehensive  discussion  of  the  circulation  features  and  data  sets.  The  main  point  is  that  the  NCOM  simulation 
employed  here  adequately  approximates  the  main  features  of  the  circulation. 

Twfo  release  locations  arc  chosen  (Fig.  7b)  away  from  coastal  regions  to  avoid  for  as  long  as  possible  the 
flow  regime  change  that  can  be  induced  by  the  strong  boundary  currents  characteristic  of  the  circulation  field. 
These  mid-basin  regions  in  the  southern  and  northern  parts  of  the  Adriatic  domain  are  square  in  shape  with 
sides  of  length  30.6  km,  in  w?hich  121  drifters  are  launched  in  a  Cartesian  array.  The  domain  size  was  also 
chosen  for  an  adequate  coverage  to  evaluate  the  decorrelation  timescales  and  the  rms  velocity  fluctuations 
of  both  the  model  and  corrected  trajectories.  Since  the  flow  is  highly  variable,  two  different  release  periods, 
October  2-17  and  October  12-27,  2002,  are  considered  in  addition  to  the  different  locations.  The  parameters 


Fig.  7,  (a)  Surface  mean  flow  and  (b)  the  eddy  kinetic  energy  distribution  and  the  two  drifter  launch  areas  in  NCOM. 
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Tabic  2 

Parameters  of  modeled  trajectories  released  in  two  different  locations  and  two  different  periods  in  NCOM 


Location:  south  middle  basin 

Location:  north  middle  basin 

Period;  October  2-17,  2002 

Set-A: 

Set-C: 

cr^  “  5.7  cm  s “l;  —  42  cm  s"3 

a“m  -  3.8  cm  s'1;  <7^  -  2.5  cm  s'1 

rJJj  =■  2. 1  days;  x'm=  1.7  days 

=  1.3  days;  Tjj,  =  1.5  days 

Period:  October  12-27,  2002 

Set-B; 

Set-D: 

=  5.3  cm  s“‘ ;  <rvm  =  4.7  cm  s_1 

a"n  ~  4.3  cm  s“ 1 ;  =  2.9  cm  s' 1 

=  1.5  days;  =  2.9  days 

1.3  days;  z*m  =  1.0  days 

of  the  model  trajectories  estimated  from  these  releases  are  summarized  in  Table  2.  First,  we  observe  that  there 
is  a  significant  anisotropy  in  the  parameters,  and  second,  they  exhibit  quite  a  bit  of  variation  in  time  and  space, 
thus  providing  a  total  of  four  test  beds,  denoted  Sct-A  to  Sct-D,  to  explore  the  performance  of  the  LSGS 
model.  The  experiments  with  the  LSGS  model  are  carried  out  for  the  main  experimental  matrix  (Table  1) 
in  each  of  the  four  test  beds  {Table  2),  on  the  basis  of  121  drifters,  each  initiated  with  10  different  t/(0)  real¬ 
izations  according  to  the  ij  variance  defined  by  Eq.  (8J  in  Section  2  (namely,  a  total  of  8  x  4  x  12  x  10  =  38, 720 
synthetic  drifters).  The  NCOM  velocity  fields  have  been  low-pass  filtered  using  a  1.5-day  window  in  order  to 
remove  high-frequency  motions,  such  as  the  inertia!  oscillations  strongly  present  in  the  interior  gyres. 

Trajectories  of  particles  released  in  Set-A  (Table  2)  are  shown  for  Exp-1  to  Exp-8  in  Fig.  8.  Certain  sym¬ 
metries  are  clearly  apparent  from  a  visual  inspection  of  the  trajectories.  There  is  not  a  significant  difference 
between  the  trajectories  obtained  in  Exp-3,  Exp-6,  and  the  model,  while  those  from  Exp-4  and  Exp-2  (Exp- 
5  and  Exp-7)  exhibit  somewhat  less  (more)  dispersion  with  respect  to  the  initial  launch  positions,  and  the 
dispersion  in  Exp-1  (Exp-8)  is  significantly  less  (more)  than  is  the  case  with  the  original  model.  Qualitatively 
similar  and  consistent  behavior  is  seen  for  the  case  of  Set-D  (Fig.  9)  and  the  other  two  sets  (Sets-B  and  C.  not 
shown).  The  preliminary  conclusion  from  these  plots  is  that  the  LSGS  model  has  a  significant  effect  on  the 
particle  trajectories  that  is  consistent  with  the  behavior  expected  from  the  modification  of  a  and  t  in  Exp-1 
to  Exp-8.  In  particular,  we  note  that  the  LSGS  model  acts  successfully  not  only  to  enhance  the  particle  dis¬ 
persion,  but  also  to  reduce  it. 

The  accuracy  of  the  LSGS  model  in  modifying  the  particle  transport  as  specified  is  explored  using  the  scat¬ 
ter  plots  of  the  targeted  (~L,^;)  and  achieved  (^,^)  parameters  for  Exp-1  to  Exp-8.  Since  there  are  four  sets 
of  experiments  in  NCOM,  parameters  (<rc,  zcf  averaged  over  all  Sets  A  to  D  are  shown.  The  agreement 
between  (ac,  tc)  and  (<rr,  Tr)  is  satisfactory  in  both  the  zonal  and  meridional  directions  and  for  all  experiments 
(Fig,  10). 

Quantitative  results  for  the  dispersion  of  the  particles  with  respect  to  the  initial  launch  location  p(t)  arc 
shown  in  Fig.  11.  Based  on  the  formula  given  in  Piterbarg  (2001a),  one  would  expect  pr(t)/pm{t) « 
0,rv'Tr/(°'m\Am)  the  so-called  inertial  regime  for  large  t.  This  corresponds  to  a  factor  of  |  for  Exp-1,  ^  for 
Exp-2  and  Exp-4,  1  for  Exp-3  and  Exp-6,  2  for  Exp-5  and  Exp-7,  and  4  for  Exp-8.  This  asymptotic  estimate 
is  subject  to  the  assumption  of  large  f,  isotropy,  and  homogeneity,  which  are  not  strictly  valid  in  the  case  of  the 
NCOM  simulations,  but  provide  a  scale  for  the  differences  in  dispersion  that  should  be  anticipated.  In  light  of 
these  estimates,  the  results  illustrated  in  Fig.  1 1  appear  to  be  reasonably  consistent  and  quite  encouraging,  in 
particular  given  the  complexity  of  the  flow  field  regulated  by  the  geometry,  forcing  fields,  and  internal 
dynamics. 

In  light  of  these  systematic  experiments,  which  appear  to  show  that  the  LSGS  model  accurately  corrects  the 
correlation  time  scales  and  turbulent  velocity  fluctuations,  further  experiments  are  conducted  with  realistic 
values  observed  in  the  Adriatic  circulation.  A  number  of  studios  have  been  conducted  to  analyze  the  large  drif¬ 
ter  data  set  collected  by  P.M.  Poulain  in  the  Adriatic  Sea  (Poulain,  1999,  2001;  Falco  et  al.,  2000;  and  Maurizi 
el  at.,  2004).  The  circulation  is  generally  divided  into  dynamically-diiTerent  regions  in  space,  such  as  the  wes¬ 
tern  Adriatic  current  and  interior  gyres,  and  different  seasons  in  time,  and  the  estimated  values  for  <rr  and  tr 
exhibit  a  broad  range.  Here,  we  will  work  with  the  bulk  estimates  based  on  Falco  et  al,  (2000)  for  the  interior 
gyres,  namely  it“  =»  =s  10  cm  s-1  and  t“  «  rj’  w  2  days.  Two  sets  arc  chosen,  one  in  the  south  and  the  other 
in  the  north,  namely  Set-A  and  Set-C  (Table  2).  Therefore,  for  Set-A,  4  =  1.75,4-  =  2.38,4  =  1.05,4  =  0.85, 
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Fig.  8.  Fifteen-day  long  trajectories  ofl21  particles  for  Exp-1  to  Exp-8  and  Set- A  in  NCOM  (only  I  >7(0)  displayed).  Red  and  blue  points 
mark  the  initial  and  final  positions*  respectively.  The  model  run  (without  LSGS)  is  depicted  in  the  middle  panel. 
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Fig,  9 ,  Same  as  Fig,  4,  but  for  Set-D, 
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Fig.  10.  Scalier  plot  of  the  targeted  and  achieved  parameters  for  Exp- 1  to  Exp-8  calculated  from  particles  released  in 

NCOM,  and  averaged  over  Sets  A  to  D  in  {a)  zonal  and  (b)  meridional  directions. 
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c  100 


Fig.  11.  Dispersion  with  respect  to  the  launch  locations,  p(/)T  for  Exp~l  to  Exp~8  in:  {a}  Sei-A;  (b)  Sct-B;  (c)  Sct-C  and  (d)  Set-D  in 
NCOM, 


and  for  Set-C,  =  2.63,  =  4.0,  =  1.53,  4-  =  1-33.  Thus,  the  correlation  time  scale  appears  to  be  quite 

accurate  in  NCOM.  whereas  the  turbulent  velocity  fluctuations  arc  somewhat  underestimated,  which  can 
be  due  to  the  many  different  factors  discussed  in  Section  1 ,  A  factor  of  two  in  either  the  velocity  fluctuation 
or  time  scale  seems  to  be  a  typical  size  of  the  error  encountered  in  realistic  numerical  models  (e.g.,  GarraiTo 
et  al..  2001).  It  should  also  be  noted  that  the  data  can  contain  errors  due  to  subsampling  in  time  along  the 
trajectories  and/or  lack  of  data  representative  of  the  exact  simulation  period.  In  any  case,  our  primary  purpose 
is  to  provide  a  feel  of  the  typical  corrections  by  the  LSGS  model  to  the  modeled  particle  trajectories. 
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Fig.  12.  Comparison  of  7-day  trajectories  from  NCOM  without  (left  panels)  and  with  (right  panels)  the  LSGS  model  for  Sct-A  (top 
panels)  and  Set-C  (bottom  panels).  The  mean  flow  is  not  included  in  Ihe  particle  advcction*  (a)  Model  -  Set- A;  (b)  LSGS  -  Set  A: 
(c)  Model  -  Set  C;  (d)  LSGS  -  Set  C 


Comparison  of  7-day  trajectories  with  and  without  the  LSGS  model  arc  shown  in  Fig.  12,  As  expected  by 
the  ratios  of  err/<rm,  the  contribution  of  the  LSGS  model  to  the  transport  of  the  particles  is  rather  significant  in 
that  they  travel  much  longer  distances,  even  though  the  initial  part  of  the  trajectories  docs  not  appear  too 
different. 

Finally,  the  mean  flow  is  added  in  order  to  advect  the  particles  with  the  full  field.  The  results  (Fig.  13)  still 
indicate  a  significant  difference  in  those  cases  with  the  LSGS  model.  Since  the  regions  are  chosen  to  coincide 
with  a  small  ratio  of  mean  and  eddy  kinetic  energy,  this  result  is  to  be  expected.  The  contribution  of  the  LSGS 
model  in  other  regions,  such  as  the  western  Adriatic  current,  is  likely  to  be  relatively  smaller. 


5.  Summary  and  conclusions 

An  accurate  calculation  of  the  Lagrangian  transport  is  important  for  a  number  of  practical  problems,  such 
as  dispersion  of  pollutants,  biological  species,  and  sediments.  Precise  techniques  to  characterize  the  Lagran¬ 
gian  pathways  have  been  developed  in  the  context  of  dynamical  system  theory,  but  they  rely  on  the  accuracy 
of  the  Eulerian  velocity  fields,  which  are  typically  derived  from  ocean  and  coastal  models.  Yet,  these  models 
not  only  contain  a  series  of  errors  due  to  uncertainties  in  the  forcing  and  boundary  conditions,  model  physics, 
and  unresolved  space  and  time  scales,  but  they  are  also  very  expensive  as  they  are  typically  run  with  the  highest 
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Fig.  1 3.  Same  as  Fig.  9  but  with  the  mean  flow  field  included  in  particle  adveeiion.  (a)  Model  -  Sct-A;  (b)  LSGS  -  Set  A;  (cj  Mode!  -  Set  C; 
(d)  LSGS  -  Set  C 


resolution  permitted  by  the  available  computational  resources.  In  particular,  an  ensemble  of  experiments  is 
desirable  to  determine  model  uncertainty,  but  remains  prohibitive,  since  the  combinative  number  and  range 
of  the  parameters  that  affect  the  Eulerian  field  can  be  very  large.  Therefore,  the  problem  arises  of  how  to 
remove  or  reduce  the  errors  that  propagate  from  the  Eulerian  field  to  the  Lagrangian  transport. 

In  this  study,  a  simple  and  practical  method  is  developed  in  order  to  reduce  the  errors  in  the  Lagrangian 
transport,  given  the  time  evolution  of  a  velocity  field.  The  main  concept  is  to  employ  the  statistical  behavior  of 
particle  trajectories  from  observational  data  sets.  The  so-called  LSGS  method  is  developed,  which  supplies  a 
correctional  velocity  vector  tj(i)  in  the  Lagrangian  transport  equation  to  reduce  the  discrepancy  of  statistical 
parameters,  namely  the  correlation  time  scale  and  turbulent  velocity  fluctuation,  between  model-generated 
(rm,  crm)  and  observed  (rr,  <rT)  trajectories,  such  that  the  corrected  particle  trajectories  show  behavior  similar 
to  the  observed  ones  (rg  ^  ir,  ac  a*  ur). 

The  LSGS  method  is  derived  theoretically  in  the  context  of  a  random-flight  model,  in  which  trajectories 
with  parameters  i  and  a  are  simulated  directly,  namely  without  the  need  of  an  underlying  Eulerian  velocity 
field.  Then,  a  space-dependent  version  of  the  LSGS  method  is  derived  using  the  so-called  Markov  velocity  field 
model,  in  which  stochastic  velocity  fields  are  generated  as  an  idealized  representation  of  two-dimensional, 
homogenous,  stationary,  turbulent  flows.  The  LSGS  model  is  tested  in  the  context  of  the  Markov  velocity  field 
model  to  investigate  its  performance  in  reproducing  tc  ^  rr  and  ac  zz  ar  using  a  series  of  experiments  with 
different  ratios  of  -1-  and  f1-.  It  is  shown  that  satisfactorily  accurate  results  are  obtained. 
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The  same  scries  of  experiments  are  then  conducted  based  on  the  hourly  output  from  a  realistic  ocean  model, 
the  Navy  Coastal  Ocean  Model,  configured  for  the  Adriatic  Sea  with  high  spatial  (1-km)  resolution.  The 
model  is  subject  to  realistic  forcing  from  (4-km  COAMPS)  winds,  tides,  and  river  runoff,  and  coupled  to  a 
coarser  resolution  global  model.  Two  regions  are  selected  with  a  high  ratio  of  eddy  to  mean  kinetic  energy 
to  maximize  the  effect  of  turbulent  fluctuations  on  the  dispersion  of  particles  (rather  than  dispersion  by  the 
mean  shear)  and  away  from  the  boundaries,  where  the  flow  parameters  change  and  become  strongly  aniso¬ 
tropic.  Despite  the  general  complexity  of  the  flow  patterns,  reasonably  good  results  are  obtained  using 
the  LSGS.  Finally,  a  set  of  experiments  including  the  mean  flow  and  realistic  ir  and  ot  for  these  regions 
demonstrate  that  the  effect  of  the  LSGS  method  on  the  particle  trajectories  seems  significant,  despite  the  real¬ 
istic  forcing  and  realistic  model  circulation.  Given  the  difficulty  of  capturing  exactly  the  turbulent  fields,  the 
use  of  such  LSGS  models  in  particle  transport  routines  in  OGCMs  appears  to  be  a  necessary  and  promising 
avenue. 

The  LSGS  method  offers  three  significant  advantages.  First,  it  can  be  applied  in  principle  to  models  at  all 
resolutions.  Since  OGCMs  use  dissipative  eddy-viscosity  Laplacian  operators,  the  effect  of  turbulent  flows  at 
the  subgrid  scale  is  suppressed.  But  given  that  oceanic  drifters  feel  forces  at  their  own  physical  scale,  the  accu¬ 
racy  of  particle  dispersion  by  turbulent  fluctuations  reduces  with  increasing  grid  spacing.  For  instance,  the 
LSGS  model  could  be  particularly  suitable  for  large-scale  OGCMs,  in  which  the  tendency  of  rm  <  zc  and 
<rm  <  a,  is  encountered  (Garraffo  et  at.,  2001).  Second,  LSGS  models  can  act  not  only  to  enhance  rn,  and 
£7tn  but  also  to  reduce  them,  if  necessary,  from  the  observations.  Finally.  LSGS  models  can  be  applied  off-line 
using  the  archived  model  output  and  there  is  no  need  to  rerun  the  OGCM.  In  addition,  the  calculation  of  i/(f) 
is  not  a  computationally  expensive  procedure.  The  primary  requirement  for  the  LSGS  method  is  that  obser¬ 
vational  drifter  data  sets  with  sufficient  coverage  are  needed  to  determine  rr  and  o>  in  the  region  of  interest. 

The  limitations  of  the  LGSG  method  also  need  to  be  clarified.  The  present  LSGS  method  corrects  turbulent 
velocity  components  only.  Thus,  corrections  are  likely  to  be  more  effective  in  regions  where  the  turbulent  com¬ 
ponent  is  significant  with  respect  to  the  mean  velocity  component.  This  is  because  of  our  assumption  that  the 
turbulent  component  of  the  ocean  flow  field  is  usually  much  more  difficult  to  reproduce  in  numerical  models 
than  the  mean  circulation.  Nevertheless,  if  the  mean  currents  are  incorrect  in  boundary  currents  and  various 
frontal  zones,  the  present  LSGS  method  is  not  likely  to  be  very  helpful.  The  exact  locations  of  so-called  hyper¬ 
bolic  points  associated  with  the  mean  field  can  also  have  a  determining  impact  on  particle  dispersion.  Thus, 
this  method  is  aimed  primarily  to  fine-tune  the  Lagrangian  transport  from  an  accurate  numerical  model  of  the 
circulation  being  investigated,  under  the  practical  assumption  that  the  exact  reproduction  of  the  turbulent  flow 
field  can  be  a  very  cumbersome  effort  because  of  many  delicate  processes  that  need  to  be  resolved  (Section  1) 
while  it  is  generally  easier  to  improve  the  computation  of  the  mean  fields. 

Several  avenues  will  be  pursued  in  future  studies.  First,  rr  and  at  do  not  usually  remain  constant  as  particles 
travel  with  the  flow  field.  Therefore,  incorporation  of  varying  ^  and  along  the  particle  trajectories  needs  to 
be  considered,  particularly  for  semi-enclosed  seas,  such  as  the  Adriatic  Sea,  that  are  characterized  by  strong 
boundary  currents.  Second,  the  LSGS  method  can  be  used  in  the  calculation  of  the  Lagrangian  structures,  for 
instance,  from  the  local  finite  scale  Lyapunov  exponent  (FSLE,  Aurell  et  al.,  1997;  Artalc  et  al„  1997). 
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Appendix  A 

A.  I.  Proof  of  LSGS  model  formula  (5) 

Consider  the  stationary  solution  tj  of  (5)  and  proof  that  the  spectral  tensor  of  vc  defined  in  (4)  coincides  with 
that  of  vr.  Since  both  have  zero  mean,  that  would  imply  statistical  equivalence  of  vc  and  vr. 
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First,  notice  the  following  expression  for  the  spectral  tensors  of  real  velocity 
Ec(a>)  =  Br(a>)  =  G>2I+ArA;-ia>(Ar-A‘t), 


(A.l) 


where  the  star  means  complex  conjugate  and  a  similar  expression  for  the  model  tensor  £,n  with  replacing  At 
and  Ar  by  Am  and  A m,  respectively.  Matrix  Br  allows  factorization 

1  hr.u  +  icu  ft  \ 


B,  =  CrC;,  C,(iw)  =  (^hr*_ 


Q 


1  /xrG  +  itu , 


and  similar  factorization  is  valid  for 

=  C  m  - 

Obviously 

Q(x)  =  ArCt{x)A;\  P  =  AtCm(x)A-1  -  AC^A;1. 

To  verify  that  vc  and  vr  have  the  same  spectrum,  it  is  sufficient  to  check  the  relation 
Q(icu)£r(ict>)0(ict>)*  —  P(ia>)£r(ia>)P(iw)‘. 

That  can  be  done  straightforwardly  by  substituting  expressions  (A.1)-{A,4)  into  the  last  relation. 
The  following  formula  for  the  spectral  density  of  tj  follows  from  (6) 

F  tm\  =  ^  +  bl  Tl"g m 

n  C2  +  CO2  7l(l  +  T^CO2) 

By  integrating  this  expression  in  Q  one  obtains 
(zmb2  -  ca>-: 


(A.2) 

(A. 3) 
(A.4) 


T  = 


c(crm  -  1) 


which  becomes  (8)  after  accounting  for  (7).  Next  by  multiplying  (6)  consequently  by  rj  and  i/m  and  averaging 
we  get  for  covariance  RmfJ  =  ufmr} 


Bmtj  “ 1 ' 


ct\l  -  abu2m 


ac 


Using  the  last  two  expressions  one  can  easily  find  the  correlation  coefficient  rcm  which  is  reduced  to  the 
formula  in  the  main  text  after  some  algebra. 


A.2.  Markov  velocity  field  model  algorithm 

The  velocity  field  with  stream-function  covariance  ( 12)  can  be  obtained  as  a  solution  of  the  equation 

f'/r  =  *{t,xty)  (A. 5) 

driven  by  a  Gaussian  white-noise  forcing  with  covariance 

,JC|  ,V|  )</»(/]  +  r,Jt|  +  x,yt  +>>)  =  8ff^/?2T"‘^(/)exp(-r2//?2). 

The  forcing  in  (A. 5)  is  taken  as 

K 

<P(t,x,y)  =  ^(k)(4(0sin  0  +  tik(t) cos $) ,  k=  (Alti2),  <f>  =  k]x  +  kiy ,  (A. 6) 

*|-i*2—* 

where  k  is  a  wave  vector,  {IK  +  l)2  is  the  total  number  of  harmonics,  £k(0t>)k(0  arc  independent  white  noises 
in  both  /  and  k,  and  finally  <r(k)  is  an  amplitude  that  in  our  experiments  was  taken  as 

<r(k)  =  A1  exp{  —sk1), 
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where  parameters  >  0  are  related  to  Grm,£,tm  by 


Thus,  we  consider  a  zero  mean  flow  with  no  spin  (co  =  0)  excited  by  a  finite  number  of  random  harmonics. 
From  (A. 6)  it  follows  that  the  velocity  field  satisfies 

-i4u  +  f(t,x)  (A.7) 

with  forcing  components  along  x-  and  y- axes,  respectively,  given  by 

ft  =  fy  =  -  ^  A,o-(k)0k(r,  at.jv),  £?k  =  cos ^  sin0. 

k€£  keif 

Introduce  stochastic  processes  Pk(/),?k(0  satisfying 

dpk{i)/dt  +  pk{t)/x  =  CkW.  d^k(r)/dr  +  qk(t)/x  =  %{t ).  (A.8) 

Then  the  components  of  the  Eulerian  velocity  field  are  given  by 

ux  =  ^  A2ff(k)(/Jk(/)  cos0  -  r/k(f)  sin  0),  uy  =  ~^2  A|ff(k)(/^(f)cos  <j>  -  qk(t)  sin  0)  (A. 9) 

UK  UK 


as  follows  from  (A.7).  For  the  Lagrangian  velocity  we  have 
di>,/d<  =  8ux/dt  +  vxdux/dx  + 
dvy/dt  =  duy/dt  +  vx8uy/dx  +  vy8uyl8y\x^t)^=m. 
Finally,  from  (A.  10)  we  get 


dvx/dt  -  vx{—[/x  +  Eu{pk..qk))  4-  dvy/dt  =  ^(-1/t  +  EZi(pk,qk))  +  vxE2i(pk,qk), 


(A.ll) 


where 


£ii  =  ^A|A;u(k)0k,  E12  =  -  ^ A;ff(k)0k,  =  ^Ajer(k)0k,  £22  =  -£11 


k£A' 


keA 


k<E* 


and 


Qk  =  -4(0s‘n</>  -  '!k(0  COS0. 

Formula  (A.  1 1)  for  the  zonal  and  meridional  components  of  the  Lagrangian  velocity  is  used  to  generate  the 
trajectories  of  particles  driven  by  the  stream  function 
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